Phase Changes in selected Lennard-Jones Xi3_nYn clusters 
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Detailed studies of the thermodynamic properties of selected binary Lennard-Jones clusters of the 
type Xi3_nYn (where n=l,2,3) are presented. The total energy, heat capacity and first derivative of 
the heat capacity as a function of temperature are calculated by using the classical and path integral 
Monte Carlo methods combined with the parallel tempering technique. A modification in the phase 
change phenomena from the presence of impurity atoms and quantum effects is investigated. 
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I. INTRODUCTION 

Clusters, as aggregates of atoms or molecules that 
range in size from two to tens of thousands of monomer 
units, can be viewed as an intermediate state of matter 
between finite and bulk. Many of the cluster properties, 
such as structural and thermodynamic for example, are 
different from the corresponding bulk properties because 
of the large number of surface species and finite size ef- 
fects. 

One thermodynamic property of clusters that 
has received much experimentaUiSiiii^iSii and 

^^„.jj.„.|-^„^ ^8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.23.24.25.26.27.28 

attention is the "phase transition". Since phase tran- 
sitions are characteristic of bulk systems, we shall 
refer to the phase transformations in clusters as phase 
changes, adopting the language introduced by Berry 
and coworkers.^ In a bulk material, the phase transition 
from solid to liquid occurs at a definite temperature and 
the heat capacity at that temperature has a sharp (delta 
function like) peaki^ In clusters, the phase change 
occurs at a range of temperatures, lying between freezing 
and melting temperatures, and the heat capacity has a 
broad peak about the transition temperature owing to 
finite size effects. Between the freezing and melting tem- 
peratures an ensemble of clusters coexists in liquid-like 
and solid-like formsiS With some notable exceptions^, 
the melting transition takes place at lower temperatures 
than the corresponding bulk^*Si 

In previous work'^-'^ (hereafter referred to as paper 1), 
we have explored the energy landscape of the selected 
binary Lennard-Jones clusters of the type Xi3_nYn and 
examined the effect of adding n Y impurity atoms on the 
structures of their Xi3_n core. In the present work, we 
extend our studies of the given systems into the thermo- 
dynamic domain. Thermodynamic properties, such as 
the total energy, heat capacity and first derivative of the 



heat capacity as a function of temperature, are calculated 
in both the classical and quantum regime. Our goals are 
to examine the effect of perturbing atoms on the thermo- 
dynamic properties of the selected systems and to under- 
stand how the complex topology of the potential energy 
surface affects phase change phenomena. In addition, we 
examine the importance of quantum contributions to the 
phase change phenomena in these systems. 

From a computational point of view, these clusters are 
also interesting because of their complex potential en- 
ergy surface. The double-funnel character of their poten- 
tial energy surface makes them a particularly challenging 
TQ^se for Monte Carlo simulationsii^iSSiS Therefore, they 
constitute a good numerical test for measuring the effi- 
ciency of Monte Carlo techniques designed to overcome 
quasiergodicity in both classical2Si2^ and quantumS^ sim- 
ulations. 

The remainder of the paper is organized as follows. In 
Sec. m we give a brief review of the methods and the 
model potential we employ to calculate the thermody- 
namic properties for the given systems. In Sec. IIIII we 
present the results including the total energy, the heat 
capacity and the first derivative of the heat capacity as 
a function of temperature. The phase change behavior 
is characterized with the help of the probability distribu- 
tion of isomers as a function of energy. Finally, in Sec. II VI 
we summarize our findings. 



II. COMPUTATIONAL METHOD 

In the present Section, we describe the computational 
details of our studies involving binary clusters of the 
form Xi3_„Y„. We have chosen to study three systems: 
X12Y1, X11Y2, and X10Y3. The choice of the systems is 
motivated by the detailed studies of their PES and con- 
struction of the associated disconnectivity graphs in the 
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companion paper 
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A. Interaction potential 

The clusters are modeled by the total potential energy 
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where VLj{rij), the pairwise Lennard- Jones potential as 
a function of the distance between atoms i and j, is 
given by 



and Vc(ri) is the confining potential 
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Finally, Rc is the radius of the confining sphere^ while 
e governs the strength of the confining potential. The 
role of the confining potential Vc{ri) is to prevent atoms 
from permanently leaving the cluster since the cluster 
in vacuum at any finite temperature is metastable with 
respect to evaporation. 

The optimal choice of the parameter Rc for the con- 
fining potential has been discussed in recent worki^S, If 
Rc is taken to be too small, the properties of the system 
become sensitive to its choice, whereas large values of Rc 
can result in problems attaining an ergodic simulation. 
The classical and quantum Monte Carlo simulations pre- 
sented here have been performed with Rc = 2axx and 
e — exx- Since the potential energy surfaces of all the 
systems studied in the present work display a double- 
funnel character,^^ their thermodynamic properties are 
calculated with Monte Carlo methods coupled with 
the parallel tempering technique^^i^^i'^'^i'^'''i'^^i'^^i'*°i'^^i'^^ de- 
vised to tackle possible ergodicity problems. 
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In Eq. the constants and aij are the energy and 
length-scale parameters for the interaction of particles i 
and j. 

For the binary clusters, we need to specify both the 
"like" (X-X, Y-Y) as well as the "mixed" (X-Y) in- 
teractions. The "mixed" Lennard- Jones parameters are 
obtained from the "like" Lennard- Jones parameters by 
usual combination rulea^ 
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B. Classical Monte Carlo simulation 

The total energy (S), the constant volume heat ca- 
pacity (Cy), and the first derivative of the heat capacity 
with respect to the temperature {d{Cv) /dT)y for clus- 
ters consisting of N particles are calculated using the 
standard expressions 
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In the present work, we have chosen to examine se- 
lected binary Lennard-Jones clusters of the type Xi3_nYn 
where the impurity atoms are less massive than their core 
(X) counterparts. Based on our previous workf^l we have 
chosen the X-atoms to be argon {exx — 119.8 K, axx — 
3.405 A, mass = 39.948) and Y-atoms to be "neon- 
like". Specifically, we have chosen the {a = cryy j uxx-, 
e — tYY / ^xx) ratios for the impurity atoms to be 
(0.8,0.6) for X12Y1 and (0.8,0.5) for both XnYa and 
X10Y3, values that produce interesting classes of poten- 
tial energy surfaces topologies. The mass of the Y-atom 
impurity has been taken to be that of of neon (20.1797) 
in all calculations. 

In Eq. Q fi and Rem are the coordinates of the ith 
atom and the center of mass of the cluster, respectively. 
The center of mass of the cluster is given by 
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kB\dT 

= -2(3{Cv) + P\V^) + 2{Vf - 3(F2) iy)l (9) 

where fc^ is Boltzmann constant, T is the temperature, 
(3 = l/ksT, V is the potential energy and angular brack- 
ets denote the canonical averages with respect to the 
Boltzmann weight exp(-/3y(x)). 



1. Parallel tempering and sampling strategy 

The parallel tempering Monte Carlo simulations are 
carried out using a total of 48 parallel streams, each run- 
ning a replica of the system at a different temperature. 
The streams are independent and uncorrelated sequences 
of random numbers that can be generated simultaneously 
on multiple processors. In this paper, we have used the 
parallel random number generator library called Scal- 
able Parallel Random Number Generator (SPRNG)iii4 
Temperatures are generated in the range from Tmin to 
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in such a way that they are distributed by geomet- 



45 



ric progressior^ 



1 < J < M, 
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We have chosen T„iin = 0.2 K, T^ax = 50 K and M = 48. 
For the range of temperatures [T„i„, Tmaa;] = [0-2, 50] K, 
a number of M = 48 streams has produced acceptance 
probabiUties for swaps larger than 40% for all streams 
and for all simulations performed. 

Explicitly, the Monte Carlo simulation is conducted as 
follows. For each stream, a random walk is carried out 
through configuration space according to the Metropo- 
lis algorithm.'*^ The basic Monte Carlo step consists of 
attempted moves of the physical coordinates associated 
with a given particle. Each attempted move is either ac- 
cepted or rejected in accord with the Metropolis prescrip- 
tion. The maximum displacements have been adjusted in 
order to ensure a 40%-60% acceptance ratio in the Monte 
Carlo moves. The maximum displacements or step sizes 
are chosen in an analogous manner to the temperatures, 
to satisfy geometric progression. In other words, the step 
size at the temperature Tj is 



where 
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The required acceptance ratio between 40% and 60% is 

achieved for s 

48. 

We define a pass as the minimal set of Monte Carlo 
attempted moves over all particles in the system. Since 
the clusters of interest are made of 13 atoms, a pass con- 
sists of 13 basic steps. We also define a block as a set of 
ten thousand passes. The size of the block is sufficiently 
large that the block averages of the estimated quantities 



are independent for all practical purposes. The simu- 
lation is divided in two phases: an equilibration phase 
that consists of 100 blocks and an accumulation phase 
that consists of 400 blocks per temperature. 

An exchange of configurations between streams at ad- 
jacent temperatures has been attempted every 25 passes 
and it has been accepted or rejected according to the 
parallel tempering logici^^i^^i'^^i^'^'i'^^i^^i'^^i'^^i"^^ A stream 
at any given temperature attempts a swap of configura- 
tions with a stream at adjacent lower and higher temper- 
ature in succession. Because of this swapping strategy, 
the streams at minimum and maximum temperatures are 
involved in swaps only every 50 passes. 

All error bars quoted in the current work correspond 
to two standard deviations. In order to avoid the 
cluttering of data the error bars have not been shown. 
The error bars are comparable to the thickness of the 
lines drawn in the various graphs. 

C. Path integral Monte Carlo simulation 

For quantum simulations of the heat capacity, we 
have employed a reweighted Wiener-Fourier path inte- 
gral (RW-WFPI) method^S*^ and recently developed en- 
ergy and heat capacity estimators that can be numeri- 
cally implemented by finite difference schemes^^iiS Since 
methodology is fully described in the cited references in 
this section we only present their salient features. 

The quantum analogs of the total energy {E) , the con- 
stant volume heat capacity (Cy), and its first derivative 
with respect to the temperature {d{Cv)/dT)v are given 
with the following expressions 



(Cv) 


.^1 






(dz\ - 




Z 




z 





(14) 



(15) 



J 



1 

ks 



d{Cv 
dT 



/3 /dZ_ 

Z ['dp, y 



z \dp^ 



P fdZ\ 



Z \df3j, 



, (16) 



r 



where Z is a partition function of the system. The parti- 
tion function is an integral of the diagonal density matrix 
over whole configuration space 



Z = I p{x; /3)dx 



(17) 



In the random series path-integral representation, the 



density matrix can be written as follows 
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where Pfp{P) is the density matrix of a free particle [for 
the zth free particle Pfp,i{l3) = (TOi/27rS^/?)-^/^ ], dP[a] is 
the probability measure, defined on the space 
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with 
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k=l 



The vectors = (ai j., a2 fc, . . . , aN j,) are path vari- 
ables or independent identically distributed (i.i.d) stan- 
dard normal variables with a.\ j. = (axi^k, o,yi,kT o,zi,k)- 
The vector x-^ = (xi, X2, . . . , xat) represents the phys- 
ical variables of the system, while a for the ith particle 
is {}t? (3 /miY^"^ where rrii is its mass, aa.^ can be written 
as 



The reader should note that ct; = ((Xxi, fj/i, cr^i) and 
<^xi = o'yi = cr^i- Afc(M) are a set of functions defined 
in the following way. Let {Afe(T)}fc>i be a set of func- 
tions defined on the interval [0, 1] that, together with the 
constant function Ao(t) = 1, make up an orthonormal 
basis in L^[0, 1]. Then, Afc(u) is defined as 



Afc(u) 



Xk{t)dt. 



In practical applications, the series in Eq. I|18|l needs 
to be replaced by a finite sum. Within the RW-WFPI 
framework, a finite dimensional approximation to the ex- 
act density matrix in Eq. H18|l is given by the expression 
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A„,fe(u) 
for 1 < k < n and 



2 sin(fc7rM) 



A„,fc(u) = /(u) sin(fc7ru) 
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for n < k < 4n, are chosen so that to maximize the rate 
of convergence, 

/9„(x;/3) p(x;/3). 



The function f{u) is defined as& 
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Therefore, a path that starts and ends in the same con- 
figuration space X (called a thermal loop) can be written 



X(u) = X + (T 
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It turns out to be useful to define several auxiliary 
quantities^ 
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X„(x, a; 13) = pfp{/3) exp h/3C/„(x, a; /3)] , (25) 
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Then it is easy to show that, for k = 1,2,3, we have 
(3^ fd^Z" 



Z \d(3'' 

/esn rfx /j^3jv dF[a]X„(x, a; /9)^i?„(x, a; /3, e) 



/j[j3„ dx /j^3„ dP[a]X„(x, a; /?) 



(27) 



The quantities above can be evaluated by Monte Carlo 
integration, where the canonical averages are car- 
ried out with respect to the Boltzmann-like weight 
exp[- E^i Et" 1 ai,kV2 - /?t/«(x, a; /?)]. The deriva- 
fives of the quantity i?„ (x, a; (3) can be expressed in terms 
of the derivatives of the quantity t/„ (x, a; (3) 
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and 
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The derivatives with respect to e are evaluated numeri- 
cally by a finite difference approximation 
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1. Parallel tempering and sampling strategy 

The sampling strategy is similar to the one employed 
in the classical Monte Carlo simulations. Here, with each 
attempted move of the physical coordinate x^ of a parti- 
cle, we attempt to move a randomly selected path vari- 
able associated with that particle. A pass is defined as 
the minimal set of Monte Carlo attempted moves over 
all 13 particles in the system. We define a block as a set 
of ten thousand passes. Exceptions are simulations with 
32 and 64 path variables, where a block contains twenty 
thousand passes. For each simulation 400 blocks have 
been utilized. 

The maximum displacements for the physical coordi- 
nates are chosen analogously to those utilized in clas- 
sical simulations [see Eas. (|12ll3|) ]. while the maximum 
displacements for the path variables are chosen in the 
following way: Afcj — As^j at temperature Tj, where 



A is a constant chosen so that the acceptance ratio for 
each randomly selected path variable is between 40% and 
60%. 

We implement a parallel tempering procedure that is 
analogous to the one utilized in classical simulations. An 
exchange of configurations between streams at adjacent 
temperatures has been attempted every 25 passes and it 
has been accepted or rejected according to the parallel 
tempering prescription. By monitoring the acceptance 
ratios for all 48 streams we have found that values have 
been larger than 40% for all simulations performed. It 
should be noted that the exchange of configurations be- 
tween streams includes both particle (physical) coordi- 
nates as well as their associated path variables. 

The value of discretization step eo needed to evaluate 
a finite difference approximation to the derivatives with 
respect to e has been set to eo=2^^^. The order of error 
for the third derivative is 0(eo) while for the first and 
the second derivative is O(eo). We would like to mention 
that the error introduced in the evaluation of the first 
derivative of the heat capacity (total energy, heat capac- 
ity), from the finite difference approximation, is at least 
ten (one thousand) times smaller than the corresponding 
statistical error for all simulations performed. 

We end this section with a comment on the conver- 
gence and the error bars involved in the determination 
of the total energy and the heat capacity. The conver- 
gence of the total energy and the heat capacity has been 
tested with respect to the number of path variables. We 
have found the results to be converged for 4n = 32 path 
variables. Further increasing of the number of path vari- 
ables to 64 has yielded the same results within the sta- 
tistical error (two standard deviations). The error bars 
have not been shown in the graphs in order to avoid the 
cluttering of data. For all but the lowest temperatures 
the error bars are comparable to the thickness of the 
lines displayed in the graphs. For reference, at the lowest 
temperature at which the path integral simulations have 
been performed, T = 4 K, the following results have been 
obtained for X10Y3: (E) = -296.47±0.13, {Cv)/NkB = 
0.52±0.20, {d{Cv)/dT)v/NkB = -0.10±0.54. For the 
same system but at the temperature 7.22 K we have ob- 
tained: {E) = -293.74 ±0.07, {Cv)/NkB = 1.10 ±0.04, 
{d{Cv)/dT)v/NkB = 0.15 ± 0.06. 



D. Quenching procedure 

We have performed minimization or quenching of the 
parallel tempering Monte Carlo (classical and quantum) 
sampled configurations by implementing the Fletcher- 
Reeves-Polak-Ribiere version of the conjugate gradi- 
ent method using the algorithm given in Numerical 
Recipesm^ The quenches allow us to interpret structural 
transformations in the clusters that are associated with 
the peaks in the heat capacity curves as well as variations 
in the slopes of the caloric and the first derivative of the 
heat capacity curves. 
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FIG. 1: Disconncctivity graph for (a) X13, (b) X12Y1, (c) X11Y2, and (d) X10Y3. The energy scale is in units of txx- The 
(cr, e) values for paucls (a d) arc (1.0,1.0), (0.8,0.6), (0.8,0.5), and (0.8,0.5), respectively. Only branches leading to the 200 
lowest-energy minima are shown. 
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III. RESULTS AND DISCUSSION 



In this section, we present the results of classical and 
quantum parallel tempering Monte Carlo simulations 
performed on each of the X12Y1, X11Y2 and X10Y3 clus- 
ters. The details about the calculations are given in Sec- 
tion HIl. 

Figure n shows disconnectivity graphs for the studied 
system together with a disconnectivity graph for the ho- 
mogeneous X13 cluster. A description of the construc- 
tion of a disconnectivity graph is given in the companion 
papeu^i and references therein. A disconnectivity graph 
is a useful tool for the visualization of the underlying 
potential energy surface of the studied systems. ^^'^^ 

While the potential energy surface of the homoge- 
neous X13 cluster is characterized by a single funneljS 
those for X12Y1, X11Y2 and X10Y3 have a double-funnel 
structure. The global minimum of each system is la- 
beled by the number 1 and the next higher-lying inher- 
ent structure is labeled by the number 2. For the X12Y1 
[X11Y2] cluster, inherent structures 1 {E = —A2A26exx) 
[{E = -39.891exx)] and 2 [E ^ -42.392exx) [{E = 
— 39.787exx)] define two basins of similar energies sep- 
arated by a large energy barrier. For the X10Y3 clus- 
ter, two basins are associated with inherent structure 1 
[E — — 38.257exx) (inherent structures 2 and 3 belong 
to the basin defined by inherent structure 1) and inherent 
structure A [E = — 37.618exx)- The lowest energy basin 
of the X12Y1 cluster contains 30 inherent structures com- 
pared with 67 inherent structures in the basin associates 
with the second lowest minimum, inherent structure 2. 
Two distinct funnels on the PES of X11Y2 cluster, one 
associated with inherent structure 1 and another with 
inherent structure 2, contain 56 and 116 inherent struc- 
tures, respectively. Therefore, their lowest energy basin 
is slightly narrower than the second basin. On the other 
hand, two funnels on the PES of X10Y3 cluster, one as- 
sociated with inherent structure 1 and another with in- 
herent structure 4, contain 98 and 85 minima, respec- 
tively. The number of inherent structures associated with 
a basin often indicates the likelihood that the system will 
relax to the minimum at the bottom of that basin. 

Figure|21displays the caloric curves, i.e. the average to- 
tal energy of the clusters, as a function of temperature. 
Solid lines represent classical, while dashed lines repre- 
sent quantum results. The caloric curves of the homoge- 
neous X13 cluster are included for comparison purposes. 
Comparison of the classical and quantum caloric curves 
indicates that quantum effects are considerable over the 
entire temperature range. A variation in the slope of the 
caloric curves is characteristic of a phase change in the 
clusters. 

The constant volume heat capacity curves of the stud- 
ied clusters as a function of temperature obtained by clas- 
sical and quantum simulations are shown in Fig. |31 and 
Fig. ^ respectively. Again, the results for the X13 clus- 
ter are shown for comparison. Replacement of n X (Ar) 
atoms in the homogeneous clusters by n Y (Ne) atoms 





FIG. 2: Classical and quantum caloric curves for the studied 
systems. Solid (dashed) lines represent classical-C (quantum- 
Q) results. (E) is given in units of Kelvin per particle. The 
number of path variables employed in the path integral sim- 
ulations is 4n=32. 




T(K) 

FIG. 3: Classical heat capacities per particle for studied 
systems in units of ks as a function of temperature. 
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FIG. 4: Quantum heat capacities per particle for studied 
systems in units of ks as a function of temperature. The 
number of path variables employed in the simulations is 4n — 
32. 



TABLE I: Classical heat capacity peak parameters for 13 
atom clusters. The temperature is given in Kelvin while the 
heat capacity is given in units of fcs per particle. The values 
are obtained by a cubic spline interpolation of the parallel 
tempering data shown in Fig. |3 The heat capacity error bars 
estimates are averages of two standard deviations of the points 
near the peaks. 





Lower temperature peak Higher temperature peak 




Tpcak 


{Cv)/kbN 




(Cv)/ktN 


Xl3 






34.09±0.18 


8.27±0.06 


X12Y1 


0.65±0.01 


3.48±0.02 


31.67±0.08 


7.54±0.05 


X11Y2 


2.21±0.07 


3.36±0.02 


28.18±0.24 


6.73±0.04 


X10Y3 


11.09±0.17 


4.03±0.02 


25.61±0.05 


6.03±0.04 



manifests itself in Fig.|31in such a way that the maxima of 
the heat capacity curves are shifted toward lower temper- 
atures with respect to the maximum of the heat capacity 
of the homogeneous cluster. A similar trend is seen for 
the quantum heat capacities (see Fig. 0J. From Fig. [SJ 
it can be seen that binary clusters exhibit lower tem- 
perature peaks, in addition to higher temperature peaks. 
Their numerical values and temperatures at which they 
occur are listed in Table Both higher and lower tem- 
perature maxima are to be discussed in more details in 
sections below. As can be seen from Fig. 0] lower tem- 
perature peaks are absent from the heat capacity curves 
obtained by path integral simulations, at least in the tem- 
perature range (from 4 K to 50 K) considered here. Nu- 
merical values of the quantum heat capacity peaks with 
associated temperatures are listed in Table HTl 



TABLE II: Quantum heat capacity peak parameters for 13 
atom clusters. The temperature is given in Kelvin while the 
heat capacity is given in units of fcs per particle. The values 
are obtained by a cubic spline interpolation of the parallel 
tempering data shown in Fig.|3 The number of path variables 
employed in the quantum calculations is 4n=32. The heat 
capacity error bars estimates are averages of two standard 
deviations of the points near the peaks. 

Lower temperature peak Higher temperature peak 
{Cv)/ktN r^eafc (GV)/fcftA~ 
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peak 



X13 
X12Y1 
X11Y2 
X10Y3 



33.68±0.21 
31.23±0.27 
27.88±0.30 
24.86±0.12 



7.89±0.04 
7.07±0.04 
6.16±0.04 
5.42±0.03 



Z 



u 

V 



H 

CP 



U 

V 



5 




FIG. 5: X12Y1. Panel (a) shows classical (solid line) and 
quantum (dashed line) heat capacities per particle in units 
of fcs as a function of temperature. Panel (b) shows first 
derivative of the heat capacity per particle in units of ks- 
The solid (dashed) line represents classical (quantum) results. 
The number of path variables employed in the path integral 
simulations is 4n — 32. 



A. X12Y1 

The heat capacity curves and their first derivatives ob- 
tained by classical and path integral Monte Carlo simu- 
lations are shown in Fig. |3Ja) and (b) , respectively. The 
solid lines are classical results. 
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FIG. 6: Distributions of the inherent structures for the X12 Yi 
cluster, generated by quenching configurations sampled from 
classical parallel tempering Monte Carlo simulations at tem- 
peratures above and below the solid « solid transition (1.7 K 
and 0.3 K) and above and below the solid ^ fiquid transition 
(40.2 K and 26.0 K). Note that the lines representing the pop- 
ulation of inherent structure 2 at energy -42.392exx, labeled 
"x20" and "x2" , have been reduced by a factor of twenty and 
two, respectively. 



FIG. 7: Distributions of the inherent structures for the X12 Yi 
cluster, generated by quenching configurations sampled from 
quantum parallel tempering Monte Carlo simulations at tem- 
peratures 4 K, 10.5 K and temperatures below (26.2 K) and 
above (40.3 K) the solid liquid transition. Note that the 
lines representing the population of inherent structure 2 at en- 
ergy -42.392exx, labeled "x20" and "x2", have been reduced 
by a factor of twenty and two, respectively. 



1. Classical simulation 

The heat capacity has a broad peak at a temperature 
of about 31.7 K and a narrow, low temperature peak 
at about 0.65 K. In order to identify the phase changes 
associated with the peaks in the heat capacity curves, 
the configurations generated by parallel tempering Monte 
Carlo simulations at a temperature somewhat below and 
above the temperatures of the peaks in the heat capacity 
curves have been quenched. The results of the quenches 
are shown in Fig. |S1 Quenching from the configurations 
sampled at 0.3 K yields only inherent structure 1 (global 
minimum-Y atom is in the interior of the cluster; center 
of the icosahedron). The quenches of the configurations 
sampled at 1.7 K produce predominantly inherent struc- 
ture 2 (Y atom is on the surface of the cluster), although 
there is still a small population (0.06) of inherent struc- 
ture 1. The low temperature peak in the classical heat 
capacity curve is associated with the structural transition 
from inherent structure 1 to inherent structure 2, that is, 
a solid <-> solid phase change between two basins. 

The high probabihty (w 0.94) of finding the X12Y1 



cluster to "dwell" in inherent structure 2, at a very low 
temperature of 1.7 K, might be surprising at first sight. 
The explanation can be found in the "shape" of its dis- 
connectivity graph, see Fig. Q^b). As we noticed ear- 
lier, by examining its disconnectivity graph, the energy 
landscape displays a double-funnel structure. The fun- 
nel associated with inherent structure 2 contains twice as 
many local minima than the funnel associated with in- 
herent structure 1. Even though funnel 1 is energetically 
more favorable, owing to the larger number of minima 
associated with funnel 2, the volume of the configura- 
tion space available to funnel 2 is likely to be larger than 
one available to funnel 1 and makes it entropicaly more 
favorable. 

At a temperature T = 26 K (below the Tpeafe) the 
system occupies mainly (« 0.93) inherent structure 2. 
However, a very small population of inherent structure 
1 and higher energy amorphous structures, is present. 
As the temperature increases, the population of higher 
energy amorphous structures rapidly increases and at 
T = 40.2 K they begin to dominate. Thus, the broad 
peak in the heat capacity at T = 31.67 K is associated 
with the solid <-!■ liquid phase change. 
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2. Path integral simulation 

Unlike the classical heat capacity, the quantum heat 
capacity curve has only one peak, centered near T = 
31.23 K. It can be seen from Fig. [S] and by comparing 
the numerical results from TablesHlandHTl the solid-liquid 
phase change in the "quantum" system occurs at slightly 
lower temperature than in the "classical" one. Quan- 
tum effects lower the transition temperature by 1.4%. 
Likewise, there is a decrease in the height of the heat 
capacity maximum. The lowering of the transition tem- 
perature and the height of the heat capacity peak from 
quantum effects are expected and have been documented 
in the literatureii2*2Si2S4^ These phenomena are the con- 
sequence of the zero-point motion and/or tunneling that 
effectively raise the energy of the inherent structures, 
lower the energy barriers between them and thus allow 
for easier isomerization. 

From the distribution of inherent structures at T = 4 K 
it can be seen that the system "likes" to stay in inherent 
structure 2 with the probability 1.0 . A normal mode 
analysis of inherent structures 1 and 2 shows that esti- 
mated zero-point energy of inherent structure 1 lies about 
22 cm""'^ higher than the corresponding zero-point energy 
associated with inherent structure 2. At the temperature 
of 10.5 K, there is still a very high population (« 0.99) 
of inherent structure 2 with almost a negligible popu- 
lation of inherent structure 1 0.01). We have com- 
pared the classical distribution of inherent structures at 
approximately the same temperatures (T = 4.09 K and 
T = 10.9 K) with its quantum counterpart. It has been 
found that the global minimum in the classical case is 
slightly more likely to be populated than in the quantum 
case. At the temperatures below (26.2 K) and above 
(40.3 K) the temperature at which the solid <-> liquid 
phase change occurs, the classical and c^uantum distribu- 
tions of inherent structures are similar 



B. X11Y2 

The heat capacity curves and their first derivatives ob- 
tained by classical and path integral Monte Carlo simu- 
lations are shown in Fig. |HIa) and (b) , respectively. The 
solid lines are classical results. 



1. Classical simulation 

As is the case with X12Y1, the heat capacity curve of 
the X11Y2 cluster has two peaks. A broad, high tem- 
perature peak occurs at 28.18 K and a smaller, narrower 
peak occurs at 2.21 K. Quenching of the configurations 
generated at T = 1.2 K gives only inherent structure 1 
(one Y atom is in the interior and the other is on the sur- 
face of the cluster). The quenches of the configurations 
generated at T = 4.2 K yield inherent structures 1, 2 and 
3 with the population probability of approximately 0.2, 



7 




T(K) 

FIG. 8: X11Y2. Panel (a) shows classical (solid line) and 
quantum (dashed line) heat capacities per particle in units 
of fc_B as a function of temperature. Panel (b) shows first 
derivative of the heat capacity per particle in units of fcs. 
The solid (dashed) line represents classical (quantum) results. 
The number of path variables employed in the path integral 
simulations is 4n = 32. 

0.77 and 0.03, respectively. Inherent structures 2 and 3 
belong to the same funnel. Thus, the low temperature 
peak in the heat capacity curve is associated with the 
structural transition between inherent structure 1 and 
inherent structure 2 and correspond to a solid solid 
phase change. 

The high probability of occupying inherent structure 
2 at a relatively low temperature can be explained in 
an analogous manner to that for the X12Y1 cluster. The 
funnel associated with inherent structure 2 [see Fig.QJb)] 
contains two times more minima than the one associated 
with inherent structure 1 and has a larger entropy. 

At a temperature T = 19.5 K, inherent structures 2, 3 
and 4 associated with funnel 2 dominate in the quench 
distribution. There is also a very small population of the 
amorphous (glassy) structures. At T = 39.5 K, which 
is on the high temperature side of the higher peak in 
the heat capacity curve, the inherent structures from the 
basin 2 still dominate, but there is an appreciable pop- 
ulation of the higher energy amorphous structures. This 
indicates that the cluster undergoes a phase change from 
a solid-like to a liquid-like form. 
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FIG. 9: Distributions of the inherent structures for the Xn Y2 
cluster, generated by quenching configurations sampled from 
classical parallel tempering Monte Carlo simulations at tem- 
peratures above and below the solid « solid transition (4.2 K 
and 1.2 K) and above and below the solid liquid transition 
(39.5 K and 19.5 K). 



FIG. 10: Distributions of the inherent structures for the 
X11Y2 cluster, generated by quenching configurations sam- 
pled from quantum parallel tempering Monte Carlo simula- 
tions at temperatures 4 K, 10.5 K and temperatures below 
(20.1 K) and above (40.3 K) the solid liquid transition. 



2. Path integral simulation 

In the temperature range examined in this work, the 
quantum heat capacity curve has only one peak at T = 
27.88 K. From Fig.|Slone can see that the quantum effects 
affect the phase change in the X11Y2 cluster by shift- 
ing the sohd-hquid transition temperature toward lower 
temperatures and decreasing the hight of the maximum 
of the heat capacity curve. The change in the transi- 
tion temperature with respect to the classical result is 
approximately 1.1%. 

The quenches from the configurations sampled by 
quantum parallel tempering Monte Carlo simulations at 
4 K yield inherent structures 2, 3 and 4 (associated with 
the basin 2) with the probability 0.95, 0.04 and 0.01, 
respectively. There is no measurable population of the 
global minimum at that temperature. Comparing the 
quantum distribution of the inherent structures at 4 K 
with the classical distribution at approximately the same 
temperature, T = 4.2 K, we see that there is a size- 
able population of the global minimum in the classical 
results. Why is that? An explanation is in the quan- 
tum effects, i.e. zero-point energy. We have performed 
a normal mode analysis of inherent structures 1 and 2 



and have found that the estimate of zero-point energy 
of inherent structure 1 lies about 30 cm~^ higher than 
the corresponding zero-point energy associated with in- 
herent structure 2. At T = 10.5 K, the quenching again 
yields inherent structures 2, 3 and 4 but with probability 
0.75, 0.21 and 0.04, respectively. At the temperatures 
below (20.1 K) and above (40.3 K) the temperature at 
which the solid <-> liquid phase change occurs, the classi- 
cal and quantum distributions of the inherent structures 
are similar 



C. X10Y3 

The heat capacity curves and their first derivatives ob- 
tained by classical and path integral Monte Carlo simu- 
lations are shown in Fig. IllT a) and (b) , respectively. The 
solid lines are classical results. 



1. Classical simulation 

The heat capacity has a broad peak at a temperature 
of about 25.6 K and a narrower, low temperature peak 
at about 11.1 K. Quenching of the configurations gener- 
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FIG. 11: X10Y3. Panel (a) shows classical (solid line) and 
quantum (dashed line) heat capacities per particle in units 
of ks as a function of temperature. Panel (b) shows first 
derivative of the heat capacity per particle in units of ks- 
The solid (dashed) line represents classical (quantum) results. 
The number of path variables employed in the path integral 
simulations is 4n = 32. 



ated at 5.1 K obtains inherent structures 1 and 2 with 
probabilities of 0.99 and 0.01, respectively. Both inher- 
ent structures 1 and 2 belong to the funnel associated 
with inherent structure 1. At the temperature 15.1 K, 
the quenches of the sampled configurations obtain two 
groups of inherent structures. One group, comprised of 
the low- lying inherent structures (1, 2, 3 and 6), belongs 
to the funnel associated with inherent structure 1 and 
other group, comprised of the low-lying inherent struc- 
tures (4, 5, 7, 8 and 9) that belongs to the funnel asso- 
ciated with inherent structure 4. The population proba- 
bility of the low-lying inherent structures at the bottom 
of the funnel associated with inherent structure 1 is 0.22 
while for those that belong to the funnel associated with 
inherent structure 4 is 0.77. There is also a very small 
population of the higher-energy, glassy-like structures. 
The low temperature heat capacity peak can be char- 
acterized as a phase change arising from the structural 
transition between two groups of the inherent structures 
rather than the structural transition between individual 
inherent structures, which is the case for the X12Y1 and 
X11Y2 clusters. Therefore, it is characterized as a solid 
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FIG. 12: Distributions of the inherent structures for the 
X10Y3 cluster, generated by quenching configurations sam- 
pled from classical parallel tempering Monte Carlo simula- 
tions at temperatures above and below the solid ^ solid tran- 
sition (15.1 K and 5.1 K) and above and below the solid ^ 
liquid transition (36.1 K and 20.9 K). 



<-!■ solid phase change. As noticed earlier, the energy 
landscape of X10Y3 has two funnels with approximately 
the same number of associated inherent structures. As 
a consequence, it is equally likely that the system upon 
cooling enters either of the funnels. At the relatively low 
temperature of 5.1 K the system dwells in the funnel as- 
sociated with the global minimum, which is not the case 
for the X12Y1 and X11Y2 clusters as discussed earlier. 

At T = 20.9 K, the group of inherent structures that 
belongs to the funnel associated with inherent structure 4 
dominates in the quench distribution. However, there is 
a large number of amorphous structures that start to be 
populated. At T = 36.1 K the population of amorphous 
structures rapidly increases. The heat capacity peak at 
T — 25.61 K is, therefore, associated with the solid ^ 
liquid phase change. 



2. Path integral simulation 

Unlike the classical heat capacity, in the range of tem- 
peratures examined in the current work, the quantum 
heat capacity curve has only one peak, at the tempera- 
ture of 24.86 K. This single peak behavior can be seen 
from Fig. inl and by comparing the numerical results from 
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FIG. 13: Distributions of the inherent structures for the 
X10Y3 cluster, generated by quenching configurations sam- 
pled from quantum parallel tempering Monte Carlo simula- 
tions at temperatures 5 K, 15.3 K and temperatures below 
(21.2 K) and above (36.2 K) the solid liquid transition. 



Tables HI and UTI that the quantum effects lower the tran- 
sition temperature by 2.9% with respect to the classical 
result. 

The quenches from the configurations sampled at 5.0 K 
yield two groups of inherent structures with approxi- 
mately the same probability. Group 1 contains inherent 
structures 1 and 2, and they belong to the funnel asso- 
ciated with inherent structure 1. Group 2 is made up of 
inherent structures 4, 5, 7 and 8, and they belong to the 
funnel associated with inherent structure 4. Comparing 
the quantum distribution of the inherent structures with 
its classical counterpart at the same temperature, one 
can see that quantum effects allow for an easier isomer- 
ization between two funnels. Moreover, the quantum ef- 
fects eliminate the low temperature peak in the quantum 
heat capacity curve over the studied temperature range. 
At T = 15.3 K two groups of inherent structures, group 1 
and 2, dominate in the quench distribution. When com- 
pared to the classical distribution of inherent structures 
at approximately the same temperature, one can see that 
the population probability of group 1 is smaller than in 
the classical case. At the higher temperatures, the classi- 
cal and quantum distributions of the inherent structures 
become similar and amorphous structures dominate. As 
expected, the heat capacity peak at T = 24.86 K is asso- 
ciated with the solid <-!■ liquid phase change. 



IV. CONCLUSIONS 

In this work, we have investigated the melting behav- 
ior of the selected binary Lennard- Jones clusters of the 
form Xi3_nYn in both classical and quantum regimes. 
For the classical results, the total energy, the heat ca- 
pacity and the first derivative of the heat capacity, are 
obtained by employing classical Monte Carlo method in 
conjunction with parallel tempering, while their quantum 
counterparts are obtained by employing the path integral 
(RW-WFPI) Monte Carlo method combined with paral- 
lel tempering. Quenching of the Monte Carlo sampled 
configurations permits us to identify the structural tran- 
sitions associated with the peaks in the heat capacity 
curves. 

Classical results show that all three studied systems ex- 
hibit a low temperature peak in the heat capacity curve. 
These low temperature peaks result from a structural 
transformation between two solid-like low lying, close in 
energy, inherent structures or groups of inherent struc- 
tures. In other words, low temperature peaks are a direct 
consequence of the double-funnel structure of the under- 
lying potential energy surface of the system. The high 
temperature peaks are associated with a solid to liquid 
melting transition. 

Replacing n X (core) atoms in the homogeneous cluster 
with n lighter impurity atoms Y results in the lowering 
of the melting temperature and the hight of the peak of 
the classical heat capacity when compared to the same 
quantities in the homogeneous cluster (see Fig. O . The 
quantum heat capacity shows a similar trend (see Fig.^. 
On the other hand, low temperature peaks in the classical 
heat capacity curves are shifted toward higher tempera- 
tures as the number of atoms Y increases in the clusters. 

Quantum effects are important but relatively modest 
in the whole range of temperatures except in the low tem- 
perature regime where they are more pronounced. This is 
understandable because the dominant species in the clus- 
ters are X atoms that mimic (heavier) Ar atoms. As the 
number of the Y atoms, that mimic (lighter) Ne atoms, 
increases so does the magnitude of the quantum effect. 
In general, the magnitude of the quantum effects depends 
on the particles mass, the temperature and the relative 
strengths of the like and mixed pair interaction poten- 
tial. Quantum effects do lower the melting temperatures 
in all three studied systems. Moreover, the low tempera- 
ture peaks present in the classical heat capacities of the 
clusters disappear in the quantum case (at least in the 
temperature range considered in this work). 
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